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Introduction 


Section  1:  INTRODUCTION 


The  purpose  of  this  report  is  to  adapt  our  earlier  work  in  nonlinear  shell- 
element  methodology  (4-6,19]  to  the  specific  kinds  of  problems  encountered  in 
the  analysis  of  critically-loaded  reinforced  concrete  shell  structures.  To  this 
end,  we  develop  a  resultant-type  form  of  the  basic  equations  (of  motion)  that 
is  both  appropriate  for  the  treatment  of  thick,  layered  shells  and  leads  to  an 
efficient  finite-element  implementation  as  well.  By  eliminating  most  of  the 
work  associated  with  thickness  integration,  the  approach  should  be  particu¬ 
larly  effective  when  numerous  layers  and/or  thickness  integration  points  are 
required,  as  in  the  case  of  inelastically  responding  reinforced  concrete.  Fur¬ 
thermore,  adhering  to  a  Continuum-Based  approach,  we  continue  to  employ 
pointwise  constitutive  relations  for  accurate  through-thickness  stress  calcula¬ 
tions,  though  these  are  ultimately  presented  to  the  element  computational 
routines  in  the  form  of  resultants.  Finally,  the  implementation  of  these  new 
procedures  represents  a  straightforward  extension  of  our  current  shell-element 
software  package. 

The  organization  of  the  report  is  as  follows.  In  Section  2,  our  previ¬ 
ous  work  is  reviewed,  with  the  main  intent  of  introducing  all  of  the  perti¬ 
nent  notation.  Then  in  Section  3,  the  basic  equations,  used  as  the  starting 
point  for  finite-element  discretization  in  (19],  are  generalized  to  accomodate 
thick/layered  shells.  The  impact  of  these  revisions  on  the  finite-element  im¬ 
plementation  is  presented  in  Section  4,  where  an  outline  of  the  Global  and 
Local  solution  algorithms  is  presented  as  well.  Finally,  a  brief  summary  and 
conclusions  are  provided  in  Section  5. 
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§2.1  Background 

In  reference  [4],  a  Continuum-Based  (CB)  approach  to  nonlinear  finite-element  shell  analy¬ 
sis  is  developed.  These  developments  represent  a  generalization  of  the  “degenerated-solid” 
element,  introduced  by  Ahmad  et  al.  [l]  for  linear  analysis.  The  basic  advantage  of 
the  CB  approach  is  its  convenient  avoidance  of  complicated  (and  often  ambiguous)  shell 
equations,  while  at  the  same  time  providing  a  straightforward  vehicle  for  finite-element 
discretization  —  via  C°  displacement  interpolation.  A  second  advantage  is  the  natural  in¬ 
clusion  of  transverse-shear  deformability,  which  increases  the  thickness  regime  over  which 
the  method  is  applicable. 

The  CB  approach  presented  in  [4]  extends  the  capabilities  to  include  arbitrarily  large 
rotations,  moderately  large  strains  and  a  large  class  of  (elastic/plastic)  constitutive  mod¬ 
els.  However,  due  to  the  full  three-dimensional  character  of  the  equations  —  the  thickness 
coordinate,  2,  appears  explicitly  in  the  element  arrays  —  element  formation  can  be  quite 
expensive.  This  is  especially  true  in  the  case  of  multi-layered  inelastic  shells,  where  thick¬ 
ness  integration  may  constitute  a  substantial  portion  of  the  computing  time.  Thus,  in  [19], 
the  thickness  coordinate  is  pre-integrated  out  of  the  arrays  by  making  several  “thinness” 
assumptions,  namely: 

(Al)  Normality  of  transverse  fibers  is  assumed  in  any  fixed  configuration.  (This  does 
not  preclude  transverse  shear  deformation  on  an  incremental  basis.) 

(A2)  Curvature  effects  are  approximated  by  using  the  average  (reference  surface)  met¬ 
ric  when  computing  kinematic  quantities.  (This  is  tantamount  to  first  approxi¬ 
mation  shell  theories  [2,11]  that  neglect  terms  of  order  h/R  with  respect  to  1.) 

(A3)  Taper  of  the  shell  wall  is  assumed  to  be  gradual  enough  that:  (i)  derivatives  of 
2  along  the  surface  directions  may  be  neglected  (kinematically)  and  (ii)  the  ‘zero 
normal  stress’  hypothesis  may  be  enforced  with  respect  to  the  rtfcrenct-turface 
normal.  (This  does  not  preclude  thickness  jumps.) 

These  assumptions  leads  to  resultant-oriented  equations  (and  element  arrays)  that 
involve  only  surface  integration.  We  therefore  refer  to  the  approach  as  the  Continuum- 
Based  Resultant  (CBR)  formulation  —  to  distinguish  it  from  classical  Resultant-Based 
formulations. 
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The  CBR  approach  is  significantly  more  economical  than  the  CB  approach,  but  due  to 
the  above  assumptions  the  advantages  are  biased  towards  “thin”  shells,  with  a  degradation 
of  accuracy  anticipated  a s  h/R  (current  thickness  to  minimum  radius  of  curvature)  is 
increased. 

In  Section  3,  the  economic  advantages  of  the  CBR  approach  are  extended  to  the  “thick” 
shell  regime.  The  remainder  of  the  current  section  reviews  the  basic  CB  and  CBR  shell 
equations  and  notational  pre-requisites  to  the  sequel.  For  greater  depth  and  detail,  the 
reader  is  referred  to  [19]. 


§2.2  Continuum-Based  (CB)  Shell  Theory 


We  start  with  the  Linearized  Variational  (LV)  equations  governing  the  incremen¬ 
tal  motion  of  a  general  three-dimensional  continuum  [16,17].  Note  that  a  Lagrangian 
(material-based)  description  of  the  motion  is  being  employed  and,  unless  otherwise  spec¬ 
ified,  all  vector  and  tensor  components  are  expressed  in  a  fixed  (global)  Cartesian  basis, 
with  indices  ranging  from  1-3. 

In  operator  notation,  the  LV  equations  may  be  written 


DM{ Au)  +  D7{ Au)  =  /”*( u)  -  rn‘(u)  -  X(u) 


where  7ezt  is  the  external  force  operator,  7tnt  the  internal  force  operator,  M  the  inertial 
force  (or  mass)  operator,  D7  the  linearized  internal  force  (or  stiffness)  operator,  and  DM 
(the  linearized  mass  operator)  is  identically  equal  to  M.  The  argument,  u,  represents 
the  material  displacement  vector,  and  u,  the  corresponding  acceleration  vector  —  both 
considered  to  be  functions  of  the  material  coordinates,  Xt. 


The  incremental  displacements  and  accelerations,  Au  and  Au,  emanate  from  the 
linearization  and  are  to  be  interpreted  as  excursions  from  the  known  (current)  configuration 
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associated  with  u.  Note  that  the  entire  right-hand-side  of  (2.1)  is  thus  known,  while  the 
left-hand-side  involves  the  tangent  operators  and  the  unknown  Au.  Thus  (2.1)  constitutes 
the  basis  for  subsequent  iterative  algorithms  (Section  4). 

In  the  current  report,  we  will  concentrate  on  the  internal  force  and  stiffness  operators, 
Txnt  and  DJ ,  since  the  external  force  and  mass  operators  are  relatively  unaffected  by 
the  choice  of  shell  formulation  (i.e.,  CB  or  CBR).  Furthermore,  it  is  these  operators  that 
engender  the  bulk  of  the  cost  in  subsequent  finite-element  numerical  computations.  Their 
continuum  definitions  are  given  as  follows 


J  6*T&dV  (2.2) 

v 

+  J  V(6ui)rffV{£iUi)dV  (2.3) 
v 

* . N/  ■—  — * ’ 

Pfl»» 


where  the  integrals  are  performed  over  the  current  volume  V  of  the  body,  with  associated 
spatial  coordinates,  x<. 


and 

DI(  Au) 


J‘n'(u)  = 


f  6*rQAidV 

v _ 


In  these  expressions,  Af  is  an  incremental  strain  “vector”  corresponding  to  the  tensor 


Ac,, 


and  arranged  in  the  order 


1  /  dAu,  dAu,-  \ 

2  \  dxj  +  dxi  ) 


(2.4) 


Ac  =  [Acn,  Ac22,  AC33,  2Acj2,  2Ac23,  2AcjjJ  (2.5) 


Similarly,  £  is  the  “vector”  corresponding  to  the  Cauchy  stress  tensor,  o,  and  arranged 
in  the  order 


8.  =  l<*n,  ^aa,  ®ss»  ^ia*  <^as»  <*3i  J 


(2.6) 
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The  symbol  6  in  (2.2)  and  (2.3)  indicates  a  variational  (or  virtual)  quantity,  so  that  Sc  is 
defined  in  precisely  the  same  manner  as  Ac  (2.5),  except  with  Au,  replaced  by  S ut,  the 
virtual  displacement.  (Note:  <5u  may  also  be  viewed  as  the  weighting  function  used  to 
obtain  the  variational  equations;  hence,  we  endow  it  with  essentially  the  same  properties 
as  Au.) 

Finally,  C  is  the  “matrix"  corresponding  to  the  material  response  tensor,  C.  This 
tensor  is  associated  with  the  rate  constitutive  equations  presented  in  [7]  (see  also  §4.2). 
The  corresponding  matrix  is  arranged  such  that 


dSu j  dAufc 

dxj  Cijkl  dxi 


ScTCAc 


(2.7) 
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2.1.2  Introduction  of  Shell  Hypotheses 

Next,  the  continuum  is  specialized  to  the  form  of  a  shell,  i.eM  a  body  bounded  by  two 
closely-spaced,  generally  curved  surfaces  separated  by  a  relatively  small,  variable  thickness 
dimension  h. 

To  obtain  the  Continuum-Based  shell  equations,  we  translate  this  qualitative  descrip¬ 
tion  into  the  following  simplifying  hypotheses  (i.e.,  mathematical  constraints)  which  are 
subsequently  introduced  in  the  LV  equations  (2.1) 
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the  reference-surface  is  chosen  to  coincide  with  the  mid-surface,  z  and  h  are  related  by 
z  =  fh/2,  where  f  is  a  parameter  that  ranges  linearly  from  -1  at  the  bottom-surface  to  +1 
at  the  top-surface. 

Regarding  the  kinematic  hypotheses,  x,  are  the  coordinates  of  a  point  on  the  reference- 
surface,  ii  is  the  unit  pseudo-normal  vector  (i.e,  the  current  orientation  of  a  vector  initially 
normal  to  the  reference-surface),  and  u,  and  u,  are  the  displacements  of  the  reference- 
surface  and  psuedo-normal  (relative  to  the  reference-surface),  respectively  (see  Fig.  2). 
Thus,  (2.8)  expresses  the  hypothesis  that  normals  remain  straight  via  linearity  through  the 
thickness.  The  incremental  constraint  (2.9)  says  that  the  incremental  pseudo-normal  dis¬ 
placement  vector  Au  has  no  component  parallel  to  the  pseudo-normal  itself.  The  equation 
is  therefore  an  expression  of  the  incremental  rigidity  of  the  normal.  Note  that  any  two 
components  of  Au  in  a  plane  orthogonal  to  x  may  be  associated  with  rotation  increments, 
say  A#ip  and  A02p,  as  illustrated  in  Figure  2. 

Regarding  the  static  hypotheses,  we  have  introduced  a  local  Cartesian  coordinate 
basis  at  each  point  (£, *},£),  called  the  lamina  basis  (denoted  by  the  subscript  “f).  The 
orthogonal  lamina  coordinates,  nt  =  (x\t,X2t,x$t)  or  ( xt,yi,zi ),  are  defined  such  that 
xe  and  yi  lie  in  the  plane  tangent  to  the  the  surface  z  =  constant,  with  zt  normal  to 
that  surface  (Fig.  3).  Thus,  (2.10)  enforces  the  zero  normal  stress  constraint  in  the 
current  configuration,  while  (2.11)  enforces  the  same  constraint  on  a  rate  (or,  equivalently, 
incremental)  basis.  Note  that  for  this  latter  purpose,  an  objective  rate  is  chosen,  namely 

—  the  Truesdell  rate  [17]. 

Incorporating  the  above  constraints  in  the  continuum  variational  operators  (2.2)-(2.3) 
leads  to  the  following  CB  shell  counterparts 


.  . 

;  \ i 
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Af/  —  [  A  £yyti  2AeZyt,  2AeyMt,  2Actz^  J  (2*15) 

and 

=  l  °**t*  °vvt'  °xvii  av*c  a*xt  J  (2-16) 

and  Ce  is  obtained  by  static  condensation  of  the  *33’  component  in  C*  as  a  result  of  (2.11). 

The  lamina-subscripted  variables  are  related  to  the  corresponding  global  Cartesian 
quantities  via  the  orthogonal  transformation  matrix,  L,  defined  such  that 

=  Lv  (2.17) 


and 

te  =  LtLr  (2.18) 

where  v  is  a  generic  vector  and  t  is  a  generic  second-rank  tensor. 


Finally,  note  that  the  total  kinematic  constraint  (2.8),  expressing  the  through- 
thickness  linearity  of  the  displacement  field,  is  implicit  in  (2.12)-(2.14),  while  the  in¬ 
cremental  kinematic  constraint  (2.9)  is  more  effectively  imposed  later,  on  a  discrete  (i.e., 
finite-element  nodal)  basis. 
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§2.3  Reduction  to  Resultant  (CBR)  Equations 


The  CB  shell  equations  differ  from  conventional  shell  equations  in  that  they  involve  full 
volume  integration  (for  a  discussion  of  other  differences  refer  to  [19]).  In  [19]  (and  simi¬ 
larly  in  [15]),  the  thickness  coordinate,  £,  is  pre- integrated  out  by  making  the  “thinness*1 
assumptions  reviewed  in  Section  2.1.  Without  such  assumptions,  even  though  the  displace¬ 
ments  are  linear  in  z  (see  (2.8)),  the  strains  do  not  in  general  inherit  this  nice  feature.  The 
culprit  is  the  Jacobian  matrix,  J*,  relating  curvilinear  (£,q,z)  and  Cartesian  (xt,yt,zg) 
derivatives,  i.e., 


*(•)  _  j-rd(-) 

‘  di 


(2.19) 


where 


it  = 


d*t 

di 


(2.20) 


In  general,  is  a  full  matrix  with  at  least  a  linear  dependence  on  z\  hence,  its  inverse 
engenders  a  rational  z  dependence  in  all  laminar  derivatives  via  (2.19).  This  distributes 
the  thickness  coordinate  throughout  the  variational  operators,  (2.12)-(2.14),  and  inhibits 
any  further  simplification. 


However,  under  the  assumptions  (A1)-(A3)  given  in  Section  2.1,  J*  is  replaced  by  the 
^-independent  matrix 


it 


it  o 

0  A/2 


(2.21) 


where  the  bar  indicates  evaluation  at  the  reference  surface.  Employing  (2.20),  the  varia¬ 
tional  equations  are  reduced  to  CBR  form,  wherein  the  pertinent  operators  are  redefined 


t-  A 


L 
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7int(u) 

=  J  SeTadS 

s 

(2.22) 

Z?/— "(A  u) 

=  J  £eTDAedS 

5 

(2.23) 

DJ—m  (Au) 

=  J  6g[  S  Ag,  dS 

s 

(2.24) 

Notice  that  the  volume  integrals  have  been  replaced  by  surface  integrals  and  the 
previous  stress  and  strain  measures  replaced  by  corresponding  resultant  quantities.  These 
^-independent  quantities  are  partitioned  into  membrane,  bending  and  transverse-shear 
components  as  follows 


Ae  = 


(2.25) 


(2.26) 


D  = 


S  = 


m  q 
m'  q' 
0 


(2.27) 


(2.28) 


The  resultant  quantities  Ae,  s  and  D,  are  related  to  the  pointwise  continuum  quan 
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tities  via  a  independent  partitioning  matrix,  Z,  i.e., 


s  = 


D  = 


Z(i)  Ae($,rj) 

(2.29) 

1  z 

i 

(2.30) 

J  Zr(i)  Qi  Z(i)  dz 

(2.31) 

where 


Z  = 


I  il  0 
0  0  1 


(2.32) 


(Note:  In  [19],  the  2x2  identity  matrix  in  the  lower  right  corner  of  Z  is  replaced  with  a 
diagonal  matrix  of  shear  correction  factors.) 

In  particular,  the  individual  stress  resultants  appearing  in  (2.26)  and  (2.28)  are  defined 
as  follows 


n  = 


m  = 


q  = 


m  = 


q'  = 


T  * 

°vvr  °zvt  J 

(2.33) 

l  °zztt 

°VW  °xvt  j  *  d* 

(2.34) 

°*vt\T  di 

(2.35) 

Oyyv  &zyi  J  *  di 

(2.36) 

[°**c 

a,yt  \T  zdi 

(2.37) 

and  n,  m  and  m'  are  expanded  matrix  counterparts  of  a,  03  and  m\  respectively.  For 
example 


n  =  / 

* 


°wt  j 
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Finally,  the  individual  resultant  strain-displacement  and  gradient-displacement  rela¬ 
tions  (the  latter  required  for  DTl,om)  may  be  written 


Aem 

= 

1  d£Ut 
*•  dxt  ’ 

if 

tdAut 

,  dAvt)  ,T 
+  bxi 

(2.38) 

Ae6 

= 

i  dAu/ 

dAvt 
dyt  ’ 

id  Arii 

{~KT 

+  ^)JT 

(2.39) 

AeJ 

= 

K*jg* 

+  Av*), 

idAwt 

(-fir 

+  Au<)Jr 

(2.40) 

and 

Ag,  =  Ve(Auie),  Auit  J T  (2.41) 

where  the  laminar  derivatives  are  computed  via  the  reference  surface  Jacobian,  j/.  It 
should  be  noted  that  through  these  expressions,  the  kinematic  constraint  (2.8)  now  appears 
explicitly  in  the  LV  equations. 


i 
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§2.4  Finite-Element  Discretization 

Finite-element  matrix  equations  of  motion  are  easily  obtained  from  either  the  CB  or  CBR 
shell  equations  by  employing  an  isoparametric  [4,20]  form  of  local  approximation.  Thus, 
both  the  displacements  and  coordinates  within  an  element  subdomain  are  interpolated 
from  corresponding  nodal  quantities  as  follows 


a=l 

(2.42) 

*(£.»?)  =  X^a(£,»?)x0 

a=l 

(2.43) 

Non 

Au(£,f?)  =  Na{t,ri)  Au« 

a=l 

(2.44) 

Non 

Au(£,»7)  =  ^TN<x(Z,ri)  Aua 

0=1 

(2.45) 

where  the  subscript  “a"  ranges  over  the  number  of  element  nodes,  and  for  quadrilateral 
elements,  the  Na  are  typically  chosen  as  Lagrangian  interpolation  functions  (see  Fig.  4). 
Note  that  the  interpolation  is  performed  in  a  fixed  Cartesian  basis,  usually  lamina  or  global, 
depending  on  the  context. 

Substitution  of  (2.37)-(2.40)  into  the  LV  equations  (2.1)  yields  the  corresponding 
finite-element  matrix  equations 


MAd  +  (Kmot/  +  K*eom)Ad  = 


where  M  is  the  assembled  mass  matrix  (from  X),  F*xt  is  the  assembled  external  force 
vector  (from  ?tzt),  Kmatl  and  K?eom  are  the  assembled  material  and  geometric  stiffness 
matrices,  respectively  (from  D7m*"  and  and  F*nl  is  the  assembled  internal  force 

vector  (from  /,n<).  The  incremental  displacement  vector,  Ad,  contains  all  of  the  active 
nodal  degrees-of-freedom.  As  explained  in  Section  4,  these  degrees-of-freedom  are  usually 
expressed  in  a  shell-oriented  (rather  than  global-Cartesian)  coordinate  system. 
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The  above  arrays  are  assembled  from  the  individual  element  arrays  [20],  with  the 
linkage  established  by  the  element  displacement  vectors ,  Ad*,  that  is 

Ad*  =  {Ad*}  =  {AdA} 

where  “a”  is  an  element  node  number  and  “A"  is  the  corresponding  globed  node  number. 
At  each  shell  node,  we  have  potentially  6  degrees-of-freedom,  defined  by 


'  ■  IS) 


(2.47) 


However,  during  assembly  the  element  arrays  are  suitably  transformed  to  discretely  enforce 
the  shell  incremental  kinematic  constraint  (2.9)  as  mentioned  above.  This  reduces  the 
number  of  degrees-of-freedom  per  node  to  5  —  except  where  junctures  occur  [19J. 


In  closing  this  section,  we  give  the  specific  definitions  of  the  element  arrays  F>nt<, 
Kmat,e  and  which  result  from  substituting  (2.37)-(2.40)  into  the  CBR  shell  oper¬ 

ators  (2.22)-(2.24),  as  it  will  be  useful  to  compare  these  to  the  revised  arrays  presented  in 
Section  4.  In  terms  of  nodal  contributions,  we  have 


nnt< 

II 

"cT 

OH 

m 

Si 

(2.48) 

s 

vrmatl* 

**ab 

=  f  BaDBfcdS 

J 

(2.49) 

s 

y  gtomt 

=  [\'*1 

J  U*!1  >ttl. 

(2.50) 
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In  practice,  the  B  matrix  at  a  given  (integration)  point,  is  first  defined  in  the  lamina 
system  at  that  point  such  that 

N.n 

Ae({,^)  =  £Ba,(£,tj)Ada<  (2.52) 

a=  I 

where  the  lamina  submatrix,  Ba<)  is  related  to  the  globally  attached  version  via 


Ba  =  BajLj 


(2.53) 


The  block-diagonal  matrix,  [LJ,  is  the  expanded  lamina  transformation  for  a  nodal  DOF 
vector  (2.43),  i.e., 

IL  Ol 


[LJ 


0  L 


(2.54) 


Using  the  above  definitions,  the  CBR  version  of  Ba<  in  terms  of  the  element  shape 
functions,  Na,  is  given  as  follows 


' Na,ti 

0 

0 

1 

0 

0 

O' 

0 

Na,yt 

0 

1 

0 

0 

0 

"*,vt 

Na,ze 

0 

1 

0 

0 

0 

0 

0 

0 

1 

Na,zt 

0 

0 

0 

0 

0 

1 

0 

Na,Vt 

0 

0 

0 

0 

1 

Na,yt 

Na 

0 

0 

0 

1 

Na 

0 

0 

.  0 

0 

Na,yt 

1 

0 

Na 

0. 

(2.55) 
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Finally,  the  scalar  coefficients  appearing  in  the  CBR  geometric  stiffness  matrix  (2.46) 

are  defined  as 

9.1  =  *«(*.)*  nV£(JVk)  (2-56) 

9.1  =  V*(A Ta)*  mVt(Nb)  +  V*(iVa)‘  q Nb  (2.57) 

9,1  =  +  (2.58) 

9,1  =  VtiNaym'VtiNt,)  +  NaqftVe(Nb) 

+  Vt(Nay  q'jVfc  (2.59) 
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§3.1  Motivation 

While  the  CBR  shell  assumptions  employed  in  the  previous  section  greatly  reduce  finite- 
element  computational  costs  via  the  use  of  thickness  pre-integration,  at  the  same  time 
they  restrict  the  maximum  thickness- to-radius  regime  relative  to  the  original  CB  formula¬ 
tion.  Thus,  for  thick  shells  with  substantial  curvature  (or  curvature  change),  there  is  the 
possibility  for  a  serious  degradation  of  accuracy. 

Furthermore,  for  problems  in  which  transverse-shear  effects  become  significant,  e.g., 
in  moderately  thick  shells  with  shear-flexible  layers,  it  may  be  necessary  to  resort  to  a  more 
refined  description  of  through-thickness  kinematics  than  provided  by  the  straight-normal 
approximation  employed  in  both  CB  and  CBR  formulations. 

Hence,  our  goal  in  this  section  is  as  follows.  We  would  like  to  relax  the  maximum¬ 
thickness  limitations  engendered  in  our  previous  work,  while  at  the  same  time  retaining  the 
cost-effectiveness  of  the  thickness  pre-integration  treatment  of  multi-layered  elastic/plastic 
shells. 

§3.2  Revised  Kinematics  1:  Curvature  Effects 

The  thickness- to-radius  limitation  imposed  by  the  CBR  method  may  be  directly  attributed 
to  the  Curvature  Assumption  (A2),  wherein  the  surface  metric,  or  Jacobian,  is  held  con¬ 
stant  (at  its  reference  surface  value)  through  the  thickness.  Recall  from  §2.3  that  the 
reason  for  doing  this  is  that  it  leads  to  a  linear  through-thickness  variation  of  the  incre¬ 
mental  strains,  and  hence  to  a  simple  separation  of  direct  and  moment  stress-resultants. 
Once  such  resultants  were  computed,  the  variational  (equilibrium)  equations  were  devoid 
of  the  thickness  coordinate,  z,  and  the  formation  of  finite-element  equilibrium  arrays  was 
reduced  in  scope  from  volume  to  surface  integration. 

Upon  closer  look  at  the  Jacobian,  and  how  it  influences  the  strain-displacement  re¬ 
lations,  it  becomes  apparent  that  if  we  agree  to  retain  the  other  two  CBR  assumptions, 
i.e..  Normality  (Al)  and  Taper  (A3),  then  the  Curvature  Assumption  may  be  completely 
abandoned  without  unduly  complicating  the  equations.  The  derivation  follows. 
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3.2.1  The  Surface  Jacobian 


First,  as  observed  in  [19],  after  invoking  CBR  assumptions  (Al)  and  (A3),  the  Jacobian 
matrix  decouples  into  laminar  (in-plane)  and  transverse  partitions 


_  dxt  Jr  0 

—  “TT  — 

di  0  h/2. 

where  j*,  the  surface  Jacobian,  has  the  linear  through-thickness  variation 


Jr  =  Jr  +  *Jr 


•  _  or) 

If  f 


Jr  - 


AA  A/S 

f  t 


The  inverse  of  j*  (used  to  obtain  laminar  derivatives  from  surface-coordinate  deriva¬ 
tives)  is  thus  the  following  matrix  of  rational  polynomials  in  z. 

T-i  I  IWr  .  -7c°fT]  /„J 


where 


j r1  =  jr  +«ir'  //« 


-eo/T 

J  l 


ar\ 

-f  f 


Jr 


arj  drj 

-f  f 


.  A 

are  the  transposed  cofactor  matrices  of  j*  and  j*,  respectively,  and  the  surface  Jacobian 
determinant,  7(f),  appearing  in  the  denominator  is  the  following  quadratic  polynomial 


j(z)  =  det(j{)  =  jo+Ji*  +  >2*2 
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where 


Jo 


j\ 


h 


.dxi  dyt  dxt  dyt. 

1  d£  dr)  ~  dr)  Qi  ' 

{&Xj.dyi  _  d'xt  dyt .  .dxjdyt  dxtdyt 
d£  dr)  dr)  di'  d£  dr)  dr)  d£ 

{dTt  dyt  _  dTtdyt 
d£  drj  dr)  d£ 


(3.9) 

(3.10) 

(3.11) 


REMARK  3.1 

_ co 

Note  that  j*  ,jf  ,jo,ji»  and  ji  involve  only  the  reference-surface  quantities,  x  and  x 
and  hence  are  independent  of  z. 


REMARK  3.2 

When  the  CBR  Curvature  Assumption  (A2)  is  introduced  as  well,  equation  (3.1)  is  replaced 
by  the  simpler  approximation  (2.21)  so  that 

jr1^)  -  it1  =  ir/r/io  (3.i2) 

which  is  just  the  first  term  in  a  Taylor  expansion  of  (3.5)  about  5  =  0,  and  leads  to  a  first 
approximation  shell  theory. 

Comparing  this  with  (3.5),  we  see  that  the  price  of  discarding  the  Curvature  Assump¬ 
tion  is  an  additional  5-scaled  matrix  in  the  Jacobian  numerator  and  a  quadratic  scalar 
function  of  z  in  the  denominator. 
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3.2.2  Laminar  Displacement  Derivatives 


We  now  examine  the  effect  of  employing  the  Jacobian  matrix  (3.5),  instead  of  the 
more  restrictive  CBR  version,  for  the  computation  of  laminar  derivatives,  as  required  in 
the  CB  shell  internal  force  and  stiffness  operators  (2.12)-(2.14). 

The  in-plane  laminar  derivatives  of  displacement  are  computed  from  the  corresponding 
surface-coordinate  derivatives  via 


dAu 


'e  _ 


(P) 


%_TdAuie 

Jt  If* 


(3.13) 


where  the  ( P )  superscript  restricts  our  attention  to  planar  components,  i.e., 

ZL-1 

8xe 

?L_J 

&vt 


£1 

dx 


HE)  -  SHE) 


while  the  subscript  ‘i*  in  (3.13)  may  range  from  1-3. 


Now,  since 


dAu 


dx<P) 


dAu,.  dAu,. 

- L  +  z  - - 1 


dx 


( p ) 


dx^i 


(3.15) 


we  may  apply  (3.13)  independently  to  translational  and  rotational  components.  Employing 
the  new  definition  (3.5)  leads  to 


dAu,^  _  1 


dx[P) 


and 


dAu 


'j_  _ 


dx 


(P) 


j{*) 


1 

j(z) 


jcof  dAuj.  .  ^cofdAui, 


J< 


d$ 


+  zje 


d$ 


-eo/dAuj.  .  -co/dAuj. 
jr  +*Je  d£ 


(3.16) 


(3.17) 


Finally,  substituting  (3.16)  and  (3.17)  back  into  (3.15)  yields  the  following  useful  expression 
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for  computing  laminar  displacement  derivatives 


dx{P) 


1  f  dA  uit  _ 


dAu  i.  dAu 


’ll  =2 


mwr  ldx[p)  dx[* 


dx[P) 


;t)\ 

(3.18) 


where  the  wide  bars  and  hats  above  the  displacement  derivatives  in- 

—  a 

dicate  association  with  j4  and  j*,  respectively,  i.e. 


(3.19) 

(3.20) 

(3.21) 

(3.22) 


dAu^ 

def 

j-ydAtt,, 

dx[p) 

U  d^> 

dA  Uit 

def 

__rdA  uit 

dx[P) 

U  d{(p) 

dA5t< 

dx\P) 

def 

(h\r-Td  Aut< 

wJ‘  «{*<’ i 

dA  uit 
dx[p' 

def 

(ll'\  rTdAi,t 
\*)S‘  9 £<'’> 

i(i)  =  j{z)/jo  =  1  + 


teMt) 


(3.23) 


REMARK  3.3 

The  CBR  (curvature-restricted)  case  is  recovered  by  neglecting  j*  and  replacing  j  by  1. 
The  result  is 

3A«w  dAu  ,,  .dAu,,l 

— fpf  =  — +  * — fpf  (3.24 

dx\P)  dx[p)  dx[p) 


which  is  the  basis  for  the  definitions  given  in  §2.3. 
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3.2.3  Strain-Displacement  Relations 


Equation  (3.18)  may  be  used  to  construct  similarly  partitioned  incremental  strain' 
displacement  relations. 


First  define  the  vector  of  in-plane  incremental  strain  components  as 


AcW  = 

<  Aeyyi  > 

=  J-  [a^»  +  IAfr>  +  i2^] 

j[x)  1  J 

1  2Aezyi  ) 

where 


A_*r 


=  < 


ox  i 

ipr 


dAvi  ,  dAut 


(3.25) 


(3.26) 


Ae\ 


(P)  _ 


dAut  dAut 

axt  oxt 


dAvi  ,  dA vi 

~3yi  +  ~Syt 


(3.27) 


Af<P) 


dA  u/ 

-fix? 
dA  vt 

dAvi  +  dAui 


(3.28) 


die  By  i 

Similarly,  define  a  vector  of  transverse  shear  strain  increments  having  the  same  struc¬ 
ture  as  the  in-plane  components  except  for  an  additional  contribution  from  the  normal 
derivatives  —  d{  ■  )/d*t 

dAwt  ,  dAui 
dxt 

JAw/  ,  oiivt 

fit +  ftt 


Ac<5> 


A.„€l  + 


{dAwt  ,  -dAui/ 

f ftt  +  zf£t 

dAwt  ,  zdAw/ 

Sn+’-3m 


}*(«} 


(3.29) 
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Rewriting  (3.29)  in  the  same  form  as  (3.25),  we  obtain 

A«<5’  =  ^-Ms>+IA««s>  +  iJA«'S)+iAiiS|l 

](z)  l  J 


where 

A<P  = 

t  dAwt  \ 
dxe  \ 

\  d£Wt  I 

k  dyt  J 

i 

AfP  = 

,  dAwt  ,  dA tv,  x 

1  dAwr  ,  dA wt  l 

a£<s>  = 

,  dA  wt  v 

A*P  = 

1“:) 

3.2.4  Resultant-Oriented  Strain  Measures 
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(3.30) 


(3.31) 


(3.32) 


(3.33) 


(3.34) 


To  complete  the  kinematic  development  and  prepare  for  insertion  in  the  variational 
equations,  we  separate  out  the  ^-dependence  in  (3.25)  and  (3.30)  via  a  partitioning  matrix 
and  deal  with  an  expanded  set  of  2-independent  strain  measures.  This  is  analogous  to  the 
CBR  membrane,  bending  and  shear  partitions  presented  in  §2.3. 


Thus,  for  the  planar  (membrane/bending)  components,  we  define 

At(P)  = 


where 


<  A.P 
AfP  J 
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(3.35) 


Ae<'> 


(3.36) 
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(3.37) 


(3.38) 


(3.39) 


(3.40) 


We  intend  to  employ  (3.35)  and  (3.38)  to  pre-integrate  the  LV  equations  through  the 
thickness.  However,  this  is  postponed  until  some  additional  kinematic  effects  have  been 
incorporated. 
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§3.3  Revised  Kinematics  2:  Transverse-Shear  Effects 

In  this  section,  we  introduce  a  simple  correction  to  the  transverse-shear  strain-displacement 
relations  used  in  the  preceding  sections.  The  correction  allows  for  a  predefined,  e.g., 
parabolic,  distribution  of  AeMZi  and  through  the  thickness  as  compared  with  the 

weighted  constant  value  used  in  the  conventional  Reissner-type  theory.  Hence,  more  accu¬ 
rate  calculation  of  pointwise  transverse-shear  stresses  is  expected,  which  can  be  especially 
important  in  the  case  of  thick  inelastically  deforming  shells.  Furthermore,  such  an  ap¬ 
proach  may  eliminate  the  need  for  transverse  shear  correction  factors  (in  some  cases),  and 
hence  relieve  the  user  input  burden. 


Before  proceeding,  we  note  that  this  type  of  approach  is  not  new.  It  has  been  used 
by  others  (e.g.,  [3])  with  good  experimental  correlation,  although  mainly  in  the  context  of 
thick  homogeneous  elastic  shells.  While  the  advantages  of  the  profile  correction  are  not 
as  obvious  for  inelastic /layered  shells,  the  justification  is  that  it  is  easy  to  evaluate  and 
basically  free  of  charge. 


In  light  of  the  above  explanation,  the  earlier  definitions  of  the  transverse  shear  strain 
increments  are  replaced  by 


(3.41) 


where,  as  suggested  in  [3],  we  tentatively  employ  the  parabolic  profile 

4(z  -  zmid)2 


p(z)  =  -• 


A* 


(3.42) 


in  which  zmid  is  the  distance  from  the  shell  reference  surface  to  the  mid-surface,  i.e.,  the 
eccentricity. 


The  correction  given  by  (3.42)  is  thus  consistent  with  the  important  case  of  shear- 
stress-free  boundary  conditions  on  top  and  bottom  shell  surfaces. 

Introducing  (3.41)  in  the  partitioned  strain-displacement  relations  given  in  (3.38)  sim- 
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ply  changes  the  definition  of  (3.39)  to 


Z<5> 


(l  +  P(i)),T 
r/-\  l1 


a  l  **i  l  7(5)i] 


(3.43) 
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§3.4  Revised  Variational  Equations 

By  direct  substitution  of  the  revised  kinematic  relations  given  in  (3.25)-(3.40)  into  the  CB 
variational  equations  (2. 1),(2.12)— (2.14)  we  arrive  at  what  we  shall  call  the  CBR2  form  of 
the  equations  (for  “2nd  approximation  CBR”).  The  only  operators  affected  are  the  internal 
force  and  stiffness  (material  and  geometric).  The  revised  definitions  are  presented  here. 

In  the  following  definitions,  all  surface  integrals  apply  to  the  rtjtrtnet  surface  (just 
like  in  the  CBR  versions).  The  z  dependence  of  dS  is  not  being  neglected,  but  rather 
absorbed  into  the  integrand  via  the  identity 

dS(z)  =  j(z)dS  (3.44) 

where  dS  is  the  differential  surface  area  at  the  reference  surface  and  j  is  the  normalized 
surface  Jacobian  defined  by  (3.8)  and  (3.23).  In  the  sequel,  the  symbol  dS  is  used  instead 
of  dS  for  notational  simplicity,  that  is  the  bars  are  omitted. 

3.4.1  CBR2  Internal  Force  Operator 


where 


and 


An  = 


=  < 


•r 


,<*>  = 


4S> 

■Ss| 

,(S> 

*2 


"  / 


(1  +  P(*)) 


>  dz 


g(S) 

£(S)S2 

ff(5)7  > 


dz 


(3.45) 


(3.46) 


(3.47) 


The  (P)  and  (S)  superscripts  on  the  stress  vectors  have  the  same  meaning  as  for  the 
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strain  vectors,  i.e., 


and 


(  a 


**e 


' vvt 


v  Oxyi  , 


(3.48) 


(3.49) 


and  the  vectors,  and  s<5),  constitute  membrane/ bending  and  shear  stress-resultants 
conjugate  to  the  strain  measures  Se^  and  6e^s\  respectively. 

Note  that  j{z)  from  (3.44)  has  been  factored  into  the  stress-resultant  definitions  (3.46) 
and  (3.47). 
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3.4.2  CBR2  Material  Stiffness  Operator 


(3.50) 

where  the  in-Plane  and  transverse-Shear  constitutive-resultant  matrices  are  defined  by 

(3.51) 


D(/>) 


■Q,n 

C(P,i 

c('V-) 

[l 

c,P,i’ 

£'*V 

1  3 

M 

, sym . 

C  z*  J 

dz 


D(S) 


=  / 


(1  ±  P{*)) 
3 


r£(S) 

C(s’i 

£S)i> 

c(s,7  -I 

ttS,P 

Q,S,i* 

c(s,i7 

sym. 

C<S)i‘ 

c,s,JV 

m 

c<s’7!  J 

dz 


(3.52) 


The  constitutive  tensors,  C(P)  and  Q(S\  are  the  partitions  of  Ct  corresponding  to  the 
planar  and  transverse-shear  stress /strain  components  defined  in  (3.25)  and  (3.29),  respec¬ 
tively. 


3.4.3  CBR2  Geometric  Stiffness  Operator 


where 


and 


DT"~  = 


S  =  j  Z(0)T»(Z(c)7(i)<fc 

M 

z-  _  j-f1  "  °.| 

Hi)  loo  o  ,  J 


(3.53) 


(3.54) 


(3.55) 
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Expansion  of  (3.54)  yields  the  following  ‘initial-stress’  resultant  matrix 

I"  Sji  Sj2  S13  S14 "1 
S22  S23  S24 


S  = 

where,  for  J  and  J  =  1,  2  or  3: 

S/j  = 


S33  S34 

S44. 


(3.56) 


=  [  dz 

J  if 2) 


S [4  =  J  o{S)z‘-' 


and  by  the  zero-normal-stress  condition, 

S4)4  -  J  j(z)  at»t  dz  =  0 

M 

The  partitioned  gradient  vector,  Agt,  in  (3.53)  is  defined  as  follows 

(  Vo(Aut<)  \ 


(3.57) 


(3.58) 


(3.59) 


Ag,  = 


V2(Au,f) 


(3.60) 


where  the  subscripted  gradient  operators  are  analogous  to  the  subscripted  strain  measures, 
(3.26)-(3.28),  i.e., 

Vo(Au^)  =  (3.61) 

dxt 

dAu,,  dA  u,. 

+  (3'62) 

£Au  i. 

V2(Au,J  =  — rpf  (3.63) 

dx'e  ' 

The  CBR2  shell  operators  presented  above  are  now  in  appropriate  form  for  finite- 
element  discretization. 
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§4.1  Overview 

In  this  section,  we  discuss  the  numerical  implementation  of  the  extended  CBR,  or  CBR2, 
shell  theory  using  the  same  finite-element  framework  as  employed  in  our  previous  work 
(e.g.,  [4], [19]).  The  discussion  includes  modifications  to  the  element  strain  interpolation 
matrices  reflecting  the  kinematic  revisions  presented  in  Section  3,  corresponding  formulas 
for  computing  the  element  force  and  stiffness  arrays,  a  summary  of  specific  shell-element 
types  (e.g.,  4,  9  and  16-node)  and  an  outline  of  the  solution  procedures  used  to  obtain 
displacements  and  stresses.  Section  2.4  is  recommended  as  background  reading. 


§4.2  Revised  Element  Arrays 

4.2.1  Strain-Displacement  Matrices 

Substituting  the  C°  element  displacement  approximations,  (2.37)-(2.40)  into  the  CBR2 
incremental  strain-displacement  relations,  (3.36)  and  (3.39),  leads  to  the  following  element 
counterparts. 


Ae(P)($,f?) 


{£,  ri)  Ada 


(4.1) 


and 

N.n 

_  /  n\  _  /  e»\ 

Ae(5)(£,r?)  =  ^B^5)(£,ij)Ada 

0=1 

(4.2) 

where  B ^  and  B^5)  are  element  membrane/bending  and  transverse-shear  strain-displace¬ 
ment  submatrices  at  element  node  “a” ,  respectively;  and  Ad„  is  the  element  incremental 
displacement,  or  DOF,  vector  at  node  “a”,  resolved  into  global  components  (see  (2.43)). 


As  explained  in  Section  2.4,  the  B  matrices  are  obtained  first  with  respect  to  the  local 
lamina  basis  and  then  transformed,  that  is 
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Bis)  = 

where  \L\  was  defined  in  (2.5)  such  that 


A*;, 


=  (ijAd; 


The  definitions  of  and  in  terms  of  the  element  shape  functions  follow 
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Na 
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The  wide  bars  and  hats  in  the  above  definitions  emanate  from  (3.19)-(3.22)  where 
they  were  used  to  define  lamina  displacement  derivatives  in  terms  of  the  natural  (£,f?) 
derivatives.  Thus  _ 

fe  ■  S 
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{&}  ■  **$’{£) 


_  a 

where  the  Jacobian  matrices  and  determinant  coefficients  jo  and  jj)  were  defined  in 
section  3.2.1.  They  are  computed  within  an  element  by  using  the  coordinate  interpolation 
formulas  (2.37)-(2.38)  as  shown  in  section  4.2.5. 

4.2.2  CBR2  Element  Internal-Force  Vector 


Substitution  of  (4.1)  and  (4.2)  into  the  internal-force  operator  (3.45)  over  an  individual 
element  subdomain,  leads  to  the  corresponding  element  force  vector.  A  nodal  subvector  is 
given  by  _ 


Knt‘  =  /  (BiP>Tt^>  +  Bis)rs <s>)dS* 

Js • 


(4.10) 


where  all  quantities  are  defined,  and  it  remains  only  to  numerically  perform  the  surface 
integration. 


Substitution  of  (4.1)  and  (4.2)  into  (3.50)  yields  the  element  material-stiffness  matrix. 
A  nodal  submatrix  is  defined  as 

(4.11) 

and  typically  employs  the  same  numerical  quadrature  scheme  as  used  for  the  internal-force 
vector  (4.29). 


4-4 


4-5 


§4.2  REVISED  ELEMENT  ARRAYS 


4.2.4  CBR2  Element  Geometric-Stiffness  Matrix 

For  the  geometric  stiffness,  we  may  exploit  the  simplicity  of  the  C°  interpolation  and 
obtain  (after  some  manipulation)  a  much  more  economical  form  than  suggested  by  (3.53). 
The  result  is  the  following  element  nodal  submatrix 


where 


ygeome 

db 


ffaS1  SaS1 


(4.12) 


9a  5  =  ™ar(S„VAr6  +  S12™6) 

+^«T(Sji^  +  S22VNb)  (4.13) 

9at  =  VKT(SitVNl  +  Sl3fNb  +  Sl4Nb) 

+ VNaT{S22^Nb  +  S23Vtffc  +  S  24Nb)  (4.14) 

9ai  =  (VKTS2l  +  VNaT  S3,  +  NaS4l)VNb 

+(VN?  S22  +  V^aTS32  +  NaS42)^Nb  (4.15) 

9ai  ~  (S22VNb  +  S23V7Vfr  +  S24Nb) 

+VNaT(S32VNi  +  S33VM  +  S34iVt) 

+Na{S4iVNl  +  S49Nb)  (4.16) 
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The  shape-function  gradient  operators 
and  (4.11),  i.e., 

VNl  = 

and 

VNa  = 


4.2.5  Computation  of  Lamina  Transformations  and  Jacobians 

-  A  f  f  # 

The  lamina  transformation  matrix,  L,  and  the  Jacobian  quantities,  j*,  j*,  jo,  ji  and  j2, 
appear  (or  are  required)  in  most  of  the  revised  element  arrays.  The  following  is  a  convenient 
step-by-step  procedure  for  computing  these  quantities  at  an  element  integration  point.  As 
a  prerequisite,  it  is  assumed  that  the  current  element  nodal  coordinates,  xfl,  and  the  shape 
functions  Na,  Na,^,  and  JVa,n  are  given. 


in  the  above  definitions  correspond  to  (4.10) 


f  1 
1  »*.ve  I 

{  ^  1 

1  »*,v<  I 


(4.17) 


(4.18) 


Step  0  “Normal”  Vectors  (x,).  Before  processing  any  of  the  integration  points,  element 
nodal  psuedo-normal  vectors,  $c0,  should  be  evaluated  (for  je,  j\  and  j2).  By  the 
Normality  Assumption  (Al),  these  may  be  computed  from  the  reference  surface 
coordinates,  i.e., 


Za  —  ®{a  *  ®fjo 

(4.19) 

where 

dx,e  .  %r3Nb.e._ 

~  Ua)  Xb 

e«a 

(4.20) 

®»J<» 

yrdNb,c\- 
-  -  JL  dTf 

PS  | 

(4.21) 

and  (a  refers  to  the  (£,q)  coordinates  at  node  “a".  This  requires  the 

?"-•?!  nation 

of  the  shape-function  natural  derivatives  at  nodes  as  well  as  integration  points. 
However,  this  involves  no  additional  computation,  since  natural  derivatives  are 
the  same  for  all  elements  of  a  given  type  and  hence  may  be  pre-computed. 


The  following  steps  are  then  performed  at  each  element  integration  point: 
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Step  1  Tarn 


= 

f?<« 

II 

XI 

a 

(4.22) 

«n(0  = 

- 

(4.23) 

M0  = 

—(0 

=  L  af  ({)x< 

a=l  v 

(4.24) 

Mfl  = 

— (£) 

=  E^(0S. 

*> 

(4.25) 

Step  2  Laminar  Transformations  (  L  )  .  The  lamina-global  transformation  matrix  (see 
(4.5))  is  defined  at  a  given  interior  point  as 

* 

L  -  (4.26) 

kj 

where  eZ(,  ey<  and  e,(  are  orthogonal  unit  vectors  parallel  to  the  local  z*,  ye  and 
zi  axes,  respectively.  To  obtain  an  unbiased  lamina  system  —  with  respect  to  the 
natural  (surface)  coordinates,  £,rj  —  we  construct  the  unit  vectors,  £*,,€•  v,,e*. 
as  follows: 

(a)  4,  =  «<  x  S.,/11  •  ||  (4.27) 


(4.28) 

(4.29) 

(4.30) 

(4.31) 


(«) 

=  e{  x  en/li  •  il 

W 

ex 

=  (e«  +  en)/||  •  || 

(c) 

6b 

=  ez<  x  e* 

M 

**< 

v/2 

=  ^(ex  -  6B) 

(«) 

y/2 

=  ^-(ex  +  eB) 

Step  3 


j,  =  i?  i  =  ; 


(e*<  *  «<)  t  (e«4  •  e„) 
(ey<  •  e<)  |  (6y<  •  S„) 


(4.32) 


i  ... 
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and 


1 1 


f  a*f' 

s*r  l  = 

(«*<  ‘  h«) 

(ez<  •  h„) 

L  9i 

9i)  J 

.(*¥<  *  h«) 

(4.33) 


Step  4  Jacobian  Determinants  and  Inverses.  The  final  step  is  to  compute  the  inverses  of 
j;  and  j*  and  use  formulas  (3.9)— (3.11)  to  obtain  the  determinant  coefficients  jo,j\ 
and  j2 •  Given  these  primitive  quantities,  it  is  then  straightforward  to  compute 
the  shape  function  laminar  derivatives  (4.8)-(4.9)  and  hence  all  of  the  kinematic 
quantities  required  for  the  element  force  and  stiffness  arrays. 
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§4.3  Specific  Shell-Element  Types 

The  elements  to  be  used  in  conjunction  with  the  current  (CBR2)  shell  formulation  are 
essentially  the  same  as  those  used  in  our  earlier  work.  Several  different  kinds  of  4,  9  and 
16-node  elements  are  being  adapted  for  this  purpose  (Fig.  4).  All  have  in  common  the  use 
of  Lagrange,  C°,  interpolation  and  selective/reduced  Gaussian  surface  quadrature  [9j. 

The  B  (“B-bar”)  technique  of  building  the  selective  (component-by-component)  in¬ 
tegration  into  the  definition  of  the  B-matrix  [8,19]  is  also  used  here.  However,  since  the 
CBR2  B-matrix  has  a  different  structure  than  the  CBR  version,  selective  integration  has 
a  slightly  different  connotation.  Such  differences  are  explained  in  the  following  sample 
element  descriptions. 

4.3.1  4-Node  Elements 

The  CBR2  approach  was  developed  primarily  for  use  with  curved  elements,  since  these 
are  expected  to  be  most  effective  for  large-deformation  problems  involving  thick,  inelastic 
shells  (i.e.,  for  reinforced  concrete  applications).  Since  the  4-node  elements  are  basically 
flat ,  they  do  not  benefit  from  the  curvature-corrections  introduced  in  section  3.2  —  al¬ 
though  they  may  benefit  from  the  parabolic  shear  corrections  (§3.3).  Nevertheless,  they 
will  be  employed  for  purposes  of  comparison  with  the  higher-order  elements.  In  all  cases, 
the  current  4-node  elements  use  bilinear  shape  functions  and  selective/reduced  integration 
on  internal  force,  stiffness  and  geometric  stiffness  arrays;  refer  to  [19]  for  details. 
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4.3.2  0-Node  Elements 


The  Heterosis  element  [10],  which  combines  Lagrange  (biquadratic)  interpolation  for 
rotations  and  Serendipity  (8-node)  interpolation  for  translations,  is  currently  the  recom¬ 
mended  9-node  element.  To  enhance  convergence  without  giving  rise  to  communicable 
mechanisms,  we  use  selective/reduced  integration  on  both  transverse-shear  and  membrane 
strains.  This  corresponds  to  2  x  2  quadrature  on  the  entire  matrix,  as  well  as  on 
partitions  of  that  couple  to  translational  DOFs:  columns  1-3  of  (4.6).  On  all  other 
terms,  we  use  a  normal  (i.e.,  nearly  exact)  rule  of  3  x  3,  and,  via  the  B-bar  technique, 
extrapolate  the  reduced  partitions  to  normal  quadrature  points  before  performing  the  ac¬ 
tual  integration  loop.  Furthermore,  it  has  also  been  found  effective  to  underintegrate  the 
geometric  stiffness  matrix,  as  reported  in  [19]. 

4.3.3  16-Node  Elements 

With  bi-cubic  shape  functions,  it  appears  safe  to  use  full  4x4  quadrature  on  the  16- 
node  elements,  as  neither  locking  nor  rank  deficiency  are  formally  present  [18].  However, 
for  thin  shells  we  have  found  that  the  convergence  rate  is  significantly  improved  by  using 
selective  3x3  quadrature  in  the  same  manner  as  described  above  for  the  9-node  element. 
It  is  likely  that  similar  improvement  will  be  observed  for  thicker  shell  problems. 

REMARK  4.1 

Ideally,  we  would  prefer  uniform  reduced  integration  for  all  of  the  above  elements,  as 
this  significantly  reduces  both  element  formation  and  stress  computation  time.  However, 
this  also  increases  rank  deficiency,  and  with  it  the  chances  of  activating  spurious  modes. 
Presently,  rigorous  techniques  for  correcting  rank  deficiency  (via  stabilization  matrices, 
etc.)  are  being  actively  pursued  [12-14];  the  results  will  be  incorporated  in  the  element 
library  when  appropriate. 
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§4.4  The  Global  Algorithm 

In  this  section  we  outline  the  basic  steps  involved  in  obtaining  a  solution  to  the  finite- 
element  equations  of  motion  for  a  shell  structure.  The  structural  equations  (2.41)  represent 
the  global  assembly  of  the  individual  element  arrays  defined  in  section  4.2.  The  constitutive 
(stress)  algorithm,  which  is  embedded  within  the  global  (displacement)  algorithm,  will  be 
discussed  subsequently  in  Section  4.5. 


For  purposes  of  illustration,  we  consider  the  case  of  nonlinear  dynamics,  and  em¬ 
ploy  the  generalized  Newmark  method  to  integrate  the  ODE  system  (i.e.,  for  temporal 
discretization),  and  “true-Newton”  tangent  stiffness  updates  for  nonlinear  iteration.  The 
global  algorithmic  equations  to  be  solved  at  each  nonlinear  iteration  (t  +  1)  within  each 
time  step  (n  +  1)  are  then  [6]: 


(4.34) 


for  the  iterative  displacement  increment,  followed  by  the  corrector  formulas 


d 

v 

a 


(»+D 

Cn*  1 

(*+D 

C„+ 1 

<<+l) 

Cn*  1 


*cl.  ® 

v£.. + 


(4.35) 

(4.36) 

(4.37) 


where  v  and  a  are  approximations  of  d  and  d,  respectively;  0  and  7  are  Newmark  algo¬ 
rithmic  parameters;  At  is  the  time  step  increment,  and  the  effective  stiffness  matrix  and 
residual  force  vector  are  defined  as 


and 


K* 


1 

0Ata 


M  +  K 


(4.38) 


R*  =  F*“  -  F,n‘  -  Ma 


(4.39) 


The  subscript  “C"  appearing  in  (4.34)-(4.37)  denotes  the  computational  basis,  that  is, 
the  directions  used  to  express  the  assembled  degrees-of-freedom  at  each  node.  For  shells, 
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these  are  usually  aligned  in  some  manner  with  the  local  surface  coordinates,  since  this  is 
both  convenient  for  boundary  conditions  and  for  removal  of  unnecessary  normal  rotation 
degrees-of-freedom  as  well. 

Finally,  the  meaning  of  the  symbolic  update  operator  “  ©  ” ,  which  is  different  from 
simple  addition  only  for  rotations,  is  explained  under  Task  4  below. 

To  complete  the  specification  of  the  algorithm,  the  following  predictor  formulas  are 
used  at  the  beginning  of  each  new  time  step,  n  +  1 

1  —  2/3 

=  dn  ®  (Atvn  +  — ^-At2an)  (4.40) 

=  vn  +  (1  -  'rjAtan  (4.41) 

=  0  (4.42) 

The  solve/correct  sequence  (4.34)-(4.37)  is  repeated  iteratively  until  convergence  at  a 
given  time  step,  i.e.,  when  Ad  and  R*  norms  become  acceptably  small.  Then  the  time 
step  is  advanced  (n  ♦-  n  +  1)  by  updating  the  prescribed  external  force  Fext,  employing 
the  predictors  (4.40)-(4.42),  and  so  on. 

The  shell-element-related  aspects  of  the  above  algorithm  are  now  summarized  for  a  sin¬ 
gle  global  iteration  cycle.  Given  the  current  element  nodal  coordinates  {x„  normal 
vectors  {xa}^|,,  global/computational  transformations  for  translations  [Ta],  and  rota¬ 
tions  [  T„  ] ^’|  j ,  the  stress  and  constitutive  tensor  at  element  integration  points,  { <r* 
and  {  Ce  } ^*1  j ,  respectively,  and  corresponding  resultant  quantities;  the  following  tasks  are 
performed: 


d<°> 
Qn  +  1 

v(0) 

Vn+1 

.(°> 

an+l 
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Task  1  Element  Formation.  Form  element  arrays  K*,  F‘nt*,  Fext*  and  Me,  and  trans¬ 
form  to  the  nodal  computational  bases  via 


Ke  _  rpT  Tf  t  m 

abc  *  M.aAablb 


(4.43) 


for  matrices,  and 


U»inte  _ 

rac  xara 


(4.44) 


for  vectors,  where  the  nodal  block-diagonal  transformation  matrix,  T«,  is  defined 
by 


Ta 


Ta  0^ 
0  XaTa 


(4.45) 


The  orthogonal  transformations,  Ta  and  Ta  are  assumed  to  be  user-supplicd  at 
the  beginning  of  the  analysis,  with  the  rotational  triads,  T,  updated  as  explained 
under  Task  4.  The  defining  relations  for  T,  T  and  x  are  as  follows: 


Aua  -  TaAuac  (4.46) 

Aua  =  X«Auac  =  X.TaA#oc  (4.47) 

where  the  skew-symmetric  matrix,  xa»  relating  rotation  increments  to  increments 
in  relative-displacement  of  the  normal  (Au),  requires  only  the  components  of  x„, 


i.e.. 


Xa  = 


0  XZa  -Z2a 

-*3«  0  Xia 

Xja  ~~  *lo  0 


(4.48) 


Task  2  Assembly.  Assemble  the  element  arrays  into  the  global  arrays:  Kc,  F‘cnt,  F^x( 
and  Me;  combine  to  form  K*  and  R*. 


Task  3  Incremental  Solve.  Solve  the  global  matrix  equations  (4.34)  for  (If  the 

matrix,  K*,  has  not  been  updated,  this  involves  only  forward  reduction  and  back 
substitution;  otherwise  the  matrix  is  first  factored.) 


Task  4  Global  Displacement  Update.  Update  the  nodal  configuration  according  to  (4.35)- 
(4.37).  (Here  we  define  the  symbol  “  ©  ”.)  For  translational  components  of 
displacement,  (4.35)  simply  implies 

OML';,'1  =  +  {Au,}  (4.49) 

while  for  rotational  components,  we  instead  update  the  rotational  triads  via 

=  [*4]S,Q(A#4,)  (4.50) 
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where  A0ac  >s  the  rotation  increment  at  node  “A”  (the  global  equivalent  of 
element  node  “a”)  emanating  from  the  solution  vector,  Adc;  and  the  orthogonal 
matrix  operator,  Q  (explained  thoroughly  in  [19])  is  defined  functionally  as 


Q(A  •) 


I  + 


1  4-  (A02/4)  L 


A©  +  -A©2 


(4.51) 


Note  that  A  0  is  the  skew-symmetric  matrix  corresponding  to  the  vector  A 8  (just 
as  x  corresponds  to  x),  and  A 6  is  its  magnitude. 


Task  5  Element  Displacement  Update.  Using  the  total  displacements,  u^,  and  the  ro¬ 
tational  triads,  T*,  computed  globally  in  Task  4,  we  then  localize  to  the  element 
level  and  obtain  the  new  reference  surface  displacement  and  pseudo-normal  vec¬ 
tors  as 

=  {x.li.Vl  +  {S.}'?,1’  (4-52) 

ar»d  ,  .  , 

-  [*.]£,•’ [*.]?<«•>.  (453) 

Note  that  the  pseudo-normal  vectors,  x„ ,  obtained  in  this  way  are  used  only  to 
compute  displacement  increments,  Aua,  for  strain/stress  calculations  (Task  6). 
Alternatively,  when  constructing  the  CBR2  shell-element  kinematic  interpolation 
arrays  (e.g.,  B),  we  use  the  reference  surface  normals  as  described  in  Section  3. 


Task  6  Constitutive  Algorithm.  Here,  we  employ  the  updated  element  nodal  coordi¬ 
nates  (4.52)  and  pseudo-normal  vectors  (4.53),  and  the  previous  stresses  and 
material-dependent  historical  quantities  to  compute  the  current  stresses 

and  constitutive  matrices  ci1+1\  resolved  in  the  current  lamina  basis  at  each 
element  integration  point.  The  algorithm  is  sketched  in  Section  4.5. 


Task  7  Stress /Constitutive  Resultants  —  Thickness  Integration.  This  constitutive  post¬ 
processing  step,  required  in  both  CBR  and  CBR2  shell  formulations,  is  also  dis¬ 
cussed  in  Section  4.5.  The  resultant  quantities  are  employed  directly  in  the  for¬ 
mation  of  element  arrays  for  the  next  iterative  cycle. 


GO  TO  Task  1  (i  ♦-  t  +  1). 


4-15 


§4.5  THE  CONSTITUTIVE  ALGORITHM 


§4.5  The  Constitutive  Algorithm 

In  this  section  we  outline  the  procedure  that  has  been  implemented  for  computing  stresses 
and  the  corresponding  resultants  required  as  input  to  the  element  equilibrium  arrays. 
Specific  constitutive  models  (e.g.,  elasticity,  plasticity,  ...)  will  not  be  discussed  here, 
but  rather  a  generic  algorithm  for  integrating  a  class  of  rate  constitutive  equations.  The 
algorithm  is  both  implementationally  convenient  and  numerically  accurate  for  problems 
involving  large  deformations  [7]. 

4.5.1  Stress  Computation 

The  Truesdell  rate  constitutive  equations  (upon  which  our  linearized  variational  equa¬ 
tions  are  based)  may  be  written  as 


(4.54) 


where 


o  +  <rw  - 


WO~  C  £ 


(4.55) 


and 


C  =  C(<r)  (4.56) 

and  iij  and  u/,y  are  the  symmetric  and  skew  symmetric  parts  of  the  velocity  gradient  tensor, 
dui/dxj. 

We  integrate  (4.54)  to  obtain  o  via  the  following  incrementally-objective  algorithm. 
The  details  are  discussed  in  (19). 

Given  a  converged  stress  state  at  time  (or  load)  step  n,  and  a  displacement  increment 
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between  step  n  and  the  current  global  iteration  at  step  n  +  1,  we  compute 


(4.57) 

where 

(4.58) 

and 

A««„,  =  Q2[C1„A££J]q5' 

(4.59) 

The  subscripts  in,  rn+i  and  imid  refer  to  the  particular  lamina  basis  active  at  steps 
n,  n  +  1  and  n  +  1/2,  respectively  (i.e.,  the  time  of  lamina-resolution).  Except  for  the 
incremental  quantities,  the  step  number  also  refers  to  the  time  of  evaluation.  Furthermore, 

C  =  C  +  C  (4.60) 

is  the  Truesdell-modihed  material  tensor,  and 

A  -rn«  =  1  ( t Au»  +  aAu>  ^  u  611 

11  2\dx'?ii  +  dx™d)  ^  ‘  ’ 

is  the  incremental  midpoint  strain  tensor,  where 

Au  =  un_f  i  -  un  (4.62) 

and 

xmid  =  ^(x„  +  Xn+i)  (4.63) 

Finally,  Qj  and  Q3  are  orthogonal  (rotation)  matrices  that  are  functions  of  the  lamina 
transformations  and  the  skew-symmetric  part  of  5A Ui/dxrptd.  In  the  presence  of  small 
shear  deformations  (both  transverse  and  in-plane),  Qj  and  Q3  each  approach  the  identity 
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matrix,  I.  For  the  sake  of  brevity,  we  will  assume  this  to  be  the  case  in  the  present 
discussion.  The  more  general  cue  is  considered  in  [19]. 

To  compute  Ac£**  at  a  typical  shell-element  integration  point,  it  is  merely  necessary 
to  replace  x  by  x™*  in  the  definition  of  the  B  matrices  given  in  Section  4.2.1.  (Note  that 
for  shells,  (4.63)  applies  to  both  x  and  x.)  The  in-plane  and  transverse-shear  components 
of  the  incremental  midpoint  strain  tensor  are  then  evaluated  using 


Ac^""*  =  Z(p)  j^BiP)™*Ada 


Ac£:“  =  Z(s)  £  B<5)™*Ada 


(4.64) 


(4.65) 


respectively. 


This  completes  the  definition  of  kinematic  quantities  required  for  the  shell  constitutive 
algorithm.  Notice,  however,  that  the  normal  strain  component,  (Ac™*  )33.  is  not  kine¬ 
matically  available  from  (4.64)-(4.65).  For  this,  the  static  (zero-normal-stress)  constraint 
is  used  to  extract  the  unknown  strain  through  the  constitutive  equations. 

In  the  elastic  case,  we  may  solve  (4.57)  for  (Ac™^)33  such  that  (**„+,) 33  =  0-  yielding 


=  “  +  £(C^)3j(Ac^’*)j 


(4.66) 


where  Qt  and  Ac™*  are  6x6  and  6x1  matrix/vector  counterparts  of  the  tensor  quantities, 
respectively,  arranged  so  that  the  33  lamina  component  comes  last.  Substitution  of  (4.66) 
in  (4.59)  then  satisfies  the  static  constraint  identically. 

For  the  inelastic  case,  it  is  suggested  that  (4.66)  be  employed  in  a  sub-iterative  fashion, 
so  that  the  updated  (inelastic)  material  tensor,  C,  may  be  accounted  for  in  the  calculation. 
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4.5.2  Resultants 

Given  the  pointwise  stresses  and  constitutive  coefficients,  it  is  then  necessary  to  in¬ 
tegrate  through  the  thickness  for  the  corresponding  resultant  quantities  appearing  in  the 
element  equilibrium  arrays.  Namely,  s^)  and  are  required  for  the  internal  force  vector 
(3.45),  and  for  the  material  stiffness  matrix  (3.50),  and  S  for  the  geometric 
stiffness  matrix  (3.53). 

Numerically,  the  resultant  quantities  are  generated  via  one-dimensional  Gauss  quadra¬ 
ture  in  a  piecewise  fashion  through  an  optional  sequence  of  “material”  layers.  Due  to 
the  potential  diversity  of  material  properties  and  layer  thicknesses,  a  different  number  of 
quadrature  points  may  be  used  within  each  tayer.  The  following  thickness  integration 
formula  is  used  for  general  layered  shells 


\  *M0 

i  E  "‘'i*) 

/  /(*) di  88  £ 

i  /=. 

t=l 

where  hj  is  the  thickness  of  layer  /,  w«  is  the  Gauss  integration  weight  at  point  *  within 
layer  l  and  £,  is  the  corresponding  thickness  coordinate.  The  relationship  between  the 
Gauss  integration  coordinate  ft,  which  is  usually  given  in  the  bi-unit  interval  (— l,+l)  — 
within  each  layer  —  and  the  thickness  coordinate  which  is  measured  from  the  element 
reference-surface,  is  as  follows 

l  =  i,  +  (4.68) 

where  zt  is  the  z  coordinate  at  the  middle  of  layer  l. 

For  single-layer  shells,  at  least  2  thickness  integration  points  are  required  in  the  elastic 
case;  more  for  plasticity  (e.g.,  5-7).  For  multi-layer  shells,  it  is  often  possible  to  make  due 
with  less  points  per  layer  (i.e.,  >  1),  but  this  is  dependent  on  the  relative  layer  thicknesses 
and  material  properties.  It  is  strongly  advised  that  such  problem-dependent  decisions  be 
made  on  the  basis  of  numerical  experiment. 
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matrix,  I.  For  the  sake  of  brevity,  we  will  assume  this  to  be  the  case  in  the  present 
discussion.  The  more  general  case  is  considered  in  [19]. 

To  compute  at  a  typical  shell-element  integration  point,  it  is  merely  necessary 

to  replace  x  by  xmtd  in  the  definition  of  the  B  matrices  given  in  Section  4.2.1.  (Note  that 
for  shells,  (4.63)  applies  to  both  x  and  x.)  The  in-plane  and  transverse-shear  components 
of  the  incremental  midpoint  strain  tensor  are  then  evaluated  using 

A*Z]?id  =  Z{P)Y^B{aP)midAda  (4.64) 

0=1 

and 

=  Z(5)^Bi5)m‘rfAd„  (4.65) 

0=1 

respectively. 

This  completes  the  definition  of  kinematic  quantities  required  for  the  shell  constitutive 
algorithm.  Notice,  however,  that  the  normal  strain  component,  (Af™^)s3,  is  not  kine¬ 
matically  available  from  (4.64)-(4.65).  For  this,  the  static  (zero-normal-stress)  constraint 
is  used  to  extract  the  unknown  strain  through  the  constitutive  equations. 

In  the  elastic  case,  we  may  solve  (4.57)  for  (Ac™*  )33  such  that  (*/w+ ,  )s3  =  0.  yielding 

(4i3»  =  -  («<.)33 +  £(£,. MAC.').,  (4.66) 

J-  1 

where  Qt  and  A sT%i  are  6x6  and  6x  1  matrix/ vector  counterparts  of  the  tensor  quantities, 
respectively,  arranged  so  that  the  33  lamina  component  comes  last.  Substitution  of  (4.66) 
in  (4.59)  then  satisfies  the  static  constraint  identically. 

For  the  inelastic  case,  it  is  suggested  that  (4.66)  be  employed  in  a  sub-iterative  fashion, 
so  that  the  updated  (inelastic)  material  tensor,  C,  may  be  accounted  for  in  the  calculation. 
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An  efficient  extension  of  the  nonlinear  finite-element  procedures  presented 
in  [19]  for  thin  shells  has  been  provided  herein  for  the  thicker,  layer-model 
type  of  shell  structures  encountered  in  many  reinforced  concrete  applications. 
The  current  technique  eliminates  much  of  the  work  associated  with  through¬ 
thickness  integration  present  in  the  earlier  Continuum-Based  shell-element  for¬ 
mulation  [4].  While  a  similar,  resultant-type  simplification  was  introduced  in 
[19]  to  first  approximation  (i.e.,  h/R  <  1),  the  more  general  case  is  now  han¬ 
dled  by  including  some  additional  terms  that  while  expanding  the  element 
kinematic  arrays,  still  restrict  thickness  integration  to  stress  and  constitutive 
quantities.  This  can  mean  significant  element-formation  cost  savings  for  multi¬ 
layer/elastic-plastic  shell  problems. 

In  addition,  a  higher-order  (parabolic)  transverse-shear  profile  has  been 
introduced  in  an  attempt  to  eliminate  the  ad-hoc  shear-correction  factors  in¬ 
herited  from  linear  analysis.  This  is  expected  to  be  an  improvement  for  those 
inelastic  problems  in  which  transverse  shear  plays  a  significant  role.  Neverthe¬ 
less,  the  representation  of  transverse  shear  in  such  complex  cases  will  require 
further  investigation. 

It  remains  to  evaluate  the  current  shell-element  procedures  using  some 
realistic  problems  involving  reinforced  concrete  constitutive  models.  Revisions 
to  the  element  software  implementation  are  now  in  progress;  these  include  ex¬ 
tension  of  the  element  arrays  to  account  for  thirk-shell  behavior  (as  described 
in  Section  4)  and  provisions  for  interfacing  with  appropriate  constitutive  mod¬ 
ules. 
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